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Abstract 

The concept of the configuration temperature is extended to discrete or discontinuous 
systems through the identification of the temperature as the nontrivial root of several 
integral equations regarding the distribution of the energy change upon configuration 
perturbations. 



Temperature is traditionally computed from the kinetic energy; but more recent 
studies [1-3] showed that it can be also measured from derivatives of the potential energy 
U (q) with respect to coordinates q: the reciprocal temperature f3 = l/(k B T) (with k B and 
T being the Boltzmann constant and absolute temperature, respectively) satisfies 



where " (. . .) " denotes an ensemble average, and (3 takes the respective meaning in 

different ensembles, e.g., it is the logarithmic derivative of the density of states g(E) in 
the microcanonical ensemble [i.e., /3{E) = -^logg(is) [4, 5]], or the parameter in the 
canonical one. Eq. (1) is exact in the canonical ensemble [3,5], but generally 
approximate, with an error of order l/N in a system of N degrees of freedom [2, 3]. 

Eq. (1) allows us to compute the temperature from a set of configurations without 
the knowledge of the velocities. This configuration temperature has found applications 
in, e.g., verifying a Monte Carlo simulation [2], entropic sampling [6], constructing a 
thermostat [7], free-energy calculation [8], etc. 

However, Eq. (1) requires a continuous system with up to the second derivatives 
of the potential energy; but the computation of second derivatives is often demanding, 
and the molecular potential is sometimes discontinuous by tail truncation [2, 9, 10] or 
simply discrete [9], making the formulas inapplicable. We will show below that the 
evaluation of derivatives can be replaced by that of the potential-energy change under a 
random perturbation, giving rise to more convenient expressions of the temperature. 

Given a configuration q, a small perturbation u in coordinates changes the 
potential energy by s = U(q + u) - U(q) « VU -u + \u ■ WU ■ u . Now consider a random 




(1) 




u, whose distribution is symmetric such that the probabilities of u and - u are equal, and 



component-independent such that w ( « y = <j S y (where the overline denotes an average 
over the perturbation u, and cr is the standard deviation of any u i ); then on average, we 



have, (i) e * \xx WU u = ±Wa , and (ii) s L ~ (VI/ -u) = (VU ■ VU)a to the 
leading order. We now take the ensemble averages over q (denoted by " (. . .} ", not to be 
confused with the overline) of the two; and their ratio is 



lis) (V 2 U) 



The last step follows from Eq. (1) (in fact, we have lim2^^ j (s k+ ^ « f3 for any odd 
positive integer k). It shows an asymmetry (s)>0 (since /? is usually positive), even 



though the perturbation is symmetric; this is because a potential-energy surface usually 
on average adopts a convex shape, which favors a positive s instead of a negative one 
[see Fig. 1(a)]. 

Eq. (2) is a "literal" translation of Eq. (1) in terms of a random perturbation; but it 
no longer requires derivatives of the potential energy function U (q) , and can thus be 
applied to some discrete U (q) . Nonetheless, the perturbation still needs to be small, for 
it is derived from the first two moments of the distribution of s . The following integral 
relations use all moments of the distribution, and thus allow larger perturbations. 

First consider a simple case for the canonical ensemble. If each configuration q is 
shifted by u to q' = q + u , and the potential-energy change a = U(q') -U(q), then [11] 

(exp(-/fe)) = 1 , (3) 



because 

(exp(-/k)), = Jexp(-/fe)3M!a j q = ^ = i , 

where Z(/?) = je~^ (q) dq is the configuration partition function, and exp[ ^ (q)] is the 
normalized weight for configuration q . 

Like Eq. (2), Eq. (3) determines j3 entirely from the statistical moments of s : 

XI-o^ - ^"^")^/"' = * ' ^ ut ^ ^ oes so as an i m pli c it equation. We can, however, obtain 
an explicit expression for /3 if the distribution of s is narrow or Gaussian: 

= log(exp(-^)> /? * -fi(e) p +\P 2 (te 2 ) p , 
where As = s- (e) . The only solution, apart from /3 = , is 

p, 2(e) J (As 2 ) p . (30 

Eq. (3') differs from (2) by using As instead of s on the denominator; but since {e) 
(the perturbation-averaged (^)^)and s are 0(a 2 ) and 0(a) , respectively, the 

difference is trifle for a small a . Thus, Eq. (3) can be considered as a generalization of 
Eq. (2) for larger perturbations. 

We will extend Eq. (3) to a few other ensembles below. We first show that it 
approximately applies to the temperature of potential energy fi v (U) = log g u (U) , 

where g v (U) = j*£(£/(q) - U)dq is the density of configurations of the same potential 

energy U: 

(exph/^t/)^*!. (4) 



We assume that the perturbation u has a symmetric distribution, such that u and - u are 
equally likely. Then, there is a microscopically reversibility: for a u that carries q to q' , 
the reverse -u that carries q' back to q shares the same probability. It follows that the 
flow Q>(U,U + s) , or the number of u that reach the U + s state from the U state (which 
is defined as all configurations with the potential energy being U ), equals the reverse 
flow 0(U + e,U) [see Fig. 1(b)]: 

0(U,U + s) = 0(U + s,U). (5) 
But the ratio of ®(U,U + s) to g v (U) is just p u (s) , the distribution of energy change 
s of the perturbations started from the U state; so 

Pu (e) = d>(U,U + s)/ gu (U). 
It is reasonable to assume that p u (s) changes slowly with U such that 

Pu( £ ) x Pu-e( £ ) ■ (6) 

Thus, 



Pu (-e) _ Q(U,U-s) _ Q(U-s,U) ^ gu(U~ g) 



^[-P v (U)s], (7) 



Pu (+e) Q>(U,U + s) Q>(U,U + s) g v (U) 
where we have used Eq. (5) in the second step, and Eq. (6) in the third. Finally, 

(exp(-^ u £)) u = £ exp(-/3 u s)p u (s) ds ~ p v (s) de = \, 

[see Fig. 1(c)]. This completes the argument. It also shows that Eq. (4) holds for a 
limited (and renormalized) portion of the distribution in (-£ max ,£ max ) (thus, we can drop 
large s commonly caused by clashes in perturbed configurations in a molecular system). 



Eq. (7) is stronger than Eq. (4), for it applies to any s ; geometrically it means that 
the regions A and B in Fig. 1(b) are roughly congruent. Thus, there are a few additional 
integral relations, e.g., 

(f*min{l,exp[-A,(C/)f]}) t ,«0, (8) 

for an odd k. Since min{l, exp(-/fe)} is the Metropolis acceptance probability [10, 12, 
13], Eq. (8) ( k = 1 ) means that the energy change after a Metropolis step starting from the 
U state averages to zero if the parameter /? matches P v (U) . Similarly, at the midpoint 

Pt, = \Pv (^0 ' me function exp(-/3 h s) p u (s) is also even; so 

(^exp[-i^(C/)^])^0, (9) 

with an odd integer k. In both cases, the s k can be replaced by an arbitrary odd function 

/odd 00 • 

The above relations can be extended to the canonical/microcanonical ensemble, or 
more generally, one that is composed of constant- U ensembles. In this case, we have 

1 = J Su (U)w(U)dU , where w(U) is the normalized weight of a configuration [which is 

e'P" /Z{fi) for the canonical ensemble, or (E -U) N/2 ~ l /g(E) for the microcanonical one 
[6], with g{E) being the density of states]. The distribution p w {s) in this ensemble is 

p w {s) = ( Pu (s)) w = \ ®(U,U + e)w(U)dU , 

where (...) denotes an average in this ensemble. Eq. (7) becomes 



Pj-£) (Pu(£)w(U + s)/w(U)) r i 

— ^ ^ M-ACM). - ex P [- <fi.V))A. 



where f3 w (U) = - log w(U) , and we have approximated w(U + s)/w(U) by 
Qxp[-/3 w (U)s], and ignored its correlation with Pyis) . Thus, Eqs. (4), (8) and (9), when 
used in this ensemble, give the temperature defined as (P w (JJ)) w = (-■^■logw(E/)) , 
which is f3 in the canonical ensemble, or /3(E) = (^^-j in the microcanonical one. 

Eqs. (4) and (7) resemble the Jarzynski equality [14] and fluctuation theorems 
[15]; but they are not identities in general, for the temperature determined here matches 
the logarithmic derivative of the flow 0(U,U + s) , instead of that of g v (U) , and thus 

has an 0(1/ N) error. Besides, although the integral relations allow larger perturbations 

than Eqs. (2) and (3'), £ should still be small in the thermodynamic sense that fi u (U) 

does not change too much from U to U + s . Eq. (7) also resembles the detailed balance 
condition used in the broad-histogram simulation [16]; but it was further approximated by 
Eq. (6) to allow an integral form, which avoids direct handling of the distribution. 

We use Eqs. (2)-(4), (8) and (9) in a molecular dynamics [9, 10] (MD) or Monte 
Carlo [9, 10, 12] (MC) simulation in this way: once every few trajectory frames, we 
perturb the current configuration q by a random u, and register the resulting change 
£ = U(q + u) - U(q) in the potential energy; the formulas then estimate the temperature 
(3 from the accumulated distribution p(s) . Since no derivative of U (q) is needed, they 
apply to systems with discrete or discontinuous potentials; we may further allow discrete 
configurations q and perturbations u, e.g., in the Ising model, q is a discrete spin 
configuration, and +u means to flip a random spin. 

We now discuss some numerical results. Four e-distributions p(s) from 
independent simulations on the 108-particle Lennard- Jones (LJ) fluid at p = 0.7 and 



/? = 1.0 are shown in Fig. 2(a). Except the first simulation, which used a regular MD 

(corresponding to a microcanonical-like ensemble), all simulations used the Metropolis 
algorithm (corresponding to the canonical ensemble) to generate trajectories. To avoid 
artifacts from the potential truncation, the pair potential was switched smoothly from the 
standard form u{r) = 4(r~ 12 - r~ 6 ) at r s = 2.0 to a 7th polynomial till r c = 2.5 , at which 

the polynomial and first three derivatives vanished [17]. In the first three cases, the 
perturbation applied to a random particle, and the Cartesian coordinates are displaced by 
random numbers within (-« max ,« max ) , with w max = 0.05 (the first two cases) or 0.1 (the 

third); but in the last case, it applied to all particles in the system with « max = 0.01 . The 

number of perturbations was 10 7 in each case. A comparison of the first two cases shows 
that the distribution depended little on whether the ensemble was microcanonical-like 
(MD) or canonical (MC). A comparison of the second and third cases shows that the 
amplitude of the single-particle perturbation mainly adjusted the aspect ratio of the 
distribution, but retained the generic shape, which peaked at s = . The distribution from 
the all-particle perturbation, however, adopted a different Gaussian-like shape. Despite 
the variety, the integral relations applied to all cases, e.g., the /? given by Eq. (4) were 
1.0077, 0.9990, 1.0160, and 0.9997, respectively. 

The solution process of Eq. (3) is illustrated in Fig. 2(b): we plotted (exp(-/fe)) 
against /? , and its intersection with 1.0 gave the required temperature. Three systems, 
the LJ fluid (with both single- and all-particle perturbations), 108-particle square- well [9] 
fluid and 32 x 32 Ising model were simulated in the canonical ensemble at /? = 1 .0 using 
the Metropolis algorithm. The number of perturbations was 10 7 in each case, and the 



computed /3 were 1.0160 (LJ, single), 0.9997 (LJ, all), 1.0003 (square-well), and 1.0004 
(Ising), respectively. The first two simulations on the LJ fluid were the same as the last 
two in Fig. 2(a). The square-well fluid differed from the LJ fluid by using a discrete pair 
potential, which is infinity if the particle distance r < r a , or - 1 if r a < r < r h , or 

otherwise (here r a = 1 and r h = 1.5 ); the perturbation was a single-atom displacement 

with M max = 0.1 ; p = 0.7 ; perturbations that produced clashes, hence infinite s , were 

excluded from entering Eq. (3), for clashed configurations contribute nothing to the 
partition function. The perturbation in the Ising model was a flip of a random site. 

In Fig. 2(c), the temperatures computed from Eqs. (l)-(3), (8) and (9) (k = 1 in the 
last two cases), after averaged over 10 7 repeats, were plotted against w max , the size of the 

single-particle perturbation. All formulas worked well with a small « max , but with a large 

M max > Eq. (2) and (3') quickly lost accuracy [Eq. (3') would, however, perform much 

better for the all-particle perturbation, whose distribution p{s) was Gaussian-like], while 

the integral formulas were little affected. 

The effect of the potential truncation [2] was studied using a set of MC 
simulations (10 single-particle perturbations in each, u max = 0.05 ) with the LJ potential 

deliberately truncated at various r c . The discontinuity rendered Eq. (1) approximate; and 
as shown in Fig. 2(d), Eq. (1) was indeed affected by a small r c , while Eqs. (3), (8), and 
(9) were not. 

The temperature fi v (U) of the potential energy was computed on the 32x32 Ising 
model. To cover the entire energy range, an entropic sampling [6, 18] simulation was run 
using the exact density of states g* v (U) [19] in the sampling weight for 10 7 MC steps per 



site. As shown in Fig. 3, the temperature from Eqs. (2) and (3') behaved badly except 
around U « or fi v « , while Eqs. (4), (8) and (9) showed good agreement with the 

reference, computed as 0* u (U) = -§j log g* v (U) . This is expected, as the smallest energy 
change is 4 in this case, the condition fie « 1 is hard to meet, while s is still smaller 
than the energy fluctuation to justify the use of the integral relations. But even in the 
latter case, there was an 0(1) difference between the logarithm density of states from the 

r u * 
integral log gyiU) = ftyiU'JdU' and the exact log g v (U) (see the inset of Fig. 3), 

which means that the temperature from the integral relations has an 0(1/ N) systematic 
error [as the entire potential-energy range is 0(N) ]. 

As an application, Eq. (3) can be used to efficiently match two potential energy 
functions, so that a simplified or coarse-grained function can emulate a more complex or 
realistic one. We use the following fluid example as a proof of principle. The aim was to 
make the configurations produced by the square- well potential U SCj (q) similar to those by 

the LJ potential U u (q) under p = 0.7 , f3 = 1 .0 . To do this, we first ran a simulation 
using U sq (q) as the potential energy, and then evaluated, from configurations of the 
trajectory, the effective temperature /? LJ by Eq. (3) with the single-particle perturbation 
( M max = 0- 1 ) and U (q) = U u (q) ; we repeat the process with a different parameter r a in 
U sq (q)(r h is kept at 1.5 for simplicity) until fi u matches (3. After r a was shifted from 
1.0 to 0.9657, /? LJ changed from 1.865 to 0.999; as a result (U u /N) was changed from 

-5.289 to -4.822, which was much closer to the reference -4.887 [20], although we did 
not explicitly matched the energy [we must expect a small deviation, for the method is 



approximate: the configurations generated by U sq (q) and U LJ (q) can be entirely 

different in principle]. The procedure was economical for the energy change needed by 
Eq. (3) can be evaluated locally; in contrast, had we matched the potential energy instead 
of temperature, it would involve a global evaluation, as well as a second simulation using 
U u (q) as the potential energy to compare the averages (U u ) in the two cases. 

To sum up, temperature ft can be extracted from the distribution p(s) of the 
potential-energy change e caused by configuration perturbations, as the nontrivial root of 
Eqs. (4), (8) or (9). If the perturbation or /? is small, simpler formulas, such as Eqs. (2) 

and (3'), also apply. Finally, we can generalize the formulas for an arbitrary distribution 
p(X) of X = X(q) [17], e.g., Eq. (4) reads 

(exp[-^logp(X) 

where AX* is the change of X after a MC/MD step (note, AX* =0 for a failed MC trial), 
and (...) denotes a configuration average at a fixed X in the simulation ensemble [this 

follows from the mapping: g v (U) — > p(X) , J3 V (U) — » log p(X) ]. 

The computer code of the examples can be found in 
http://simulalgo.appspot.com/rpt. 



Figure Captions 



Figure 1. (a) The intersection of the sphere (dotted) of perturbations u (isotropic, fixed 
length) with the constant- U surface at U + s has a larger circumference than that with the 
constant- U surface at U - s ; thus + s occurs more often than - s , even for a symmetric 
u. (b) The forward and reverse flows (arrows), e.g., ®(U,U + s) and <t>({7 + s,U), are 
balanced between any two macroscopic constant-t/ ensembles (ovals); and the flow 
<£>(U,U + s) is proportional to the ensemble population g v (U) . (c) The distribution 

Puis) (solid line) of the energy-change s upon perturbations. With the correct /?, 
exp(-/3s)p u (s) (dashed) coincides with the horizontal reflection p u (—s) of p v (s) , and 
thus satisfies (exp(-/fe)) = 1 ; it is not so with another ft' (dot-dashed). 

Figure 2. (a) The distributions of the potential-energy change s . (b) The intersection of 
(exp(-/fc)) an d 1.0 gives the correct /? . (c) The temperature versus the perturbation 
amplitude w max . (d) The temperature versus the potential truncation r c . 

Figure 3. The profile of the temperature of potential energy PyiU) along U in the 

f3 u (U')dU' . 
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